#ifndef ORTHAN2_HPP
#define ORTHAN2_HPP

#include "Density.hpp"
#include <armadillo>
#include <iostream>
#include "utils.hpp"

using namespace arma;
using namespace std;

namespace Orthant2{

class Orthant
{
	public:
		Orthant(mat A, mat a){ 
			_A=A;
			_a=a;
			_p=A.n_rows;
			d= new boost::math::normal_distribution<>(0,1); 
			_gamman=0;
			_i=0;
  		}
		~Orthant()
		{
			delete d;
		}
		inline double Likelihood(mat theta1)
		{
			double sum=0;
			for(int i=0;i<_p;i++)
			{
				sum+=log(pdf(*d,theta1(i,0)));

			}
			return sum;
		}
		int Set(mat eta,int i)
		{
			int viol=0;
			double a=0;
			if( i==0)
			{
				a=(_a(0,i))/_A(i,i);

			}else{

				a=(_a(0,i)-as_scalar(_A(i,span(0,i-1))*eta(span(0,i-1),0)))/_A(i,i);
			}
			//	double b=(_b(i,0)-as_scalar(_A(i,span(0,i-1))*eta(span(0,i-1),0)))/_A(i,i);
			if(eta(i,0) > a )//& (eta(i,0)<b))
			{
				viol=0;
			}else{
				viol=1;
			}
			return viol;
		}

		double Density_utn(mat theta)
		{
			return Likelihood(theta)+log(1-Set(theta,_gamman));
		}
		void Weight(mat theta,double *W)
		{
			int M=theta.n_cols;
		/*	double sum=0;
			for(int i=0;i<M;i++){
				sum+=exp(W[i]);	
			}
			for(int i=0;i<M;i++){
				W[i]=W[i]-log(sum);
		
			}

		*/	
			for(int i=0;i<M;i++)
			{
				if(exp(W[i])>0)
				{
					//cout << W[i] << " " ;
					if(Set(theta.col(i),_i)==0)
					{
						W[i]=0;
					}else{
						W[i]=-numeric_limits<double>::infinity();
					}
				}else{
					W[i]=-numeric_limits<double>::infinity();
				}		//W[i]=_Temp(i,0);
				

			}
		}
	
		inline void Set_i(int i){
			_i=i;
		}
		inline int Get_i(void){return _i;}
		inline int Get_p(void){return _p;}
		boost::math::normal_distribution<>* Get_d(void){return d;}
	private:
		mat _A;
		mat _a;
		int _i;
		//mat _b;
		double _gamman;
		int _p;
		mat _Temp;
		boost::math::normal_distribution<> *d;

};

} // end namespace
#endif
